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THE HALO MASS FUNCTION FROM EXCURSION SET THEORY. 
II. THE DIFFUSING BARRIER 
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ABSTRACT 

In excursion set theory the computation of the halo mass function is mapped into a first-passage time process 
in the presence of a barrier, which in the spherical collapse model is a constant and in the ellipsoidal collapse 
model is a fixed function of the variance of the smoothed density field. However, A^-body simulations show 
that dark matter halos grow through a mixture of smooth accretion, violent encounters and fragmentations, 
and modeling halo collapse as spherical, or even as ellipsoidal, is a significant oversimplification. In addition, 
the very definition of what is a dark matter halo, both in A^-body simulations and observationally, is a diffi- 
cult problem. We propose that some of the physical complications inherent to a realistic description of halo 
formation can be included in the excursion set theory framework, at least at an effective level, by taking into 
account that the critical value for collapse is not a fixed constant S c , as in the spherical collapse model, nor a 
fixed function of the variance a of the smoothed density field, as in the ellipsoidal collapse model, but rather is 
itself a stochastic variable, whose scatter reflects a number of complicated aspects of the underlying dynamics. 
Solving the first-passage time problem in the presence of a diffusing barrier we find that the exponential factor 
in the Press-Schechter mass function changes from exp{-(5 2 /2er 2 } to exp{-a(5 2 /2cr 2 }, where a = 1/(1 +Dg) 
and Db is the diffusion coefficient of the barrier. The numerical value of Db, and therefore the corresponding 
value of a, depends among other things on the algorithm used for identifying halos. We discuss the physical 
origin of the stochasticity of the barrier and, from recent A^-body simulations that studied the properties of the 
collapse barrier, we deduce a value Db — 0.25. Our model then predicts a ~ 0.80, in excellent agreement with 
the exponential fall off of the mass function found in A^-body simulations, for the same halo definition. Com- 
bining this result with the non-markovian corrections computed in paper I of this series, we derive an analytic 
expression for the halo mass function for gaussian fluctuations and we compare it with A^-body simulations. 

Subject headings: cosmology:theory — dark matter:halos — large scale structure of the universe 



1. INTRODUCTION 

The relation between the linear density perturbations 
at early time and the abundance of virialized dark mat- 
ter halos at the present epoch is an extremely relevant 
one in modern cosmology. In particular, primordial non- 
gaussianities leave an imprint on the abundance and on 
the clustering properties of the most massive objects, 
such as galaxy clusters, which form out of rar e fluc- 
tuations (iMatarrese et al.l Il986t iGrinstein & Wisd 119861, 
Lucchin et al.1 Il988t iMoscardini et all 1199 U iKoy a ma et al.1 
19991: IMatarrese et al.l 120001: iRobinson & Baked l200oT 
Robinson et al.1 120001) . These observational signatures are 



potentially detectable by various planned large-scale galaxy 
surveys. 

From the theoretical side, the challenge is to compute the 
number density of dark matter halos of mass M, n{M), in 
terms of the statistical properties of the primordial density 
field. The formation and evolution of dark matter halos is 
a highly complicated process, and its full dynamical com- 
plexity can only be studied by A^-body simulations. As re- 
vealed by A^-body simulations, halos grow through a messy 
mixture of violent encounters, s mooth accretion and fragmen- 
tation (see Springel et al. (2005) and the related movies at 
http://www.mpa-garching.mpg.de/galform/millennium/). 

Still, some analytic understanding of halo formation is 
highly desirable, both for obtaining a better physical intu- 
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ition, and for the flexibility under changes of models or pa- 
rameters (such as cosmological model, shape of the non- 
Gaussianities, etc.) which is the advantage of analytical 
results over very timing consuming numerical simulations. 
Presently the best available analytical technique is b ased on 
Press-Schecther (PS) theory dPress & Schechterll 19741) and its 
extension known as excursion set theory dBondet al.1 1 19911) 
and is able to reproduce, at least qualitatively, several proper- 
ties of dark matter halos seen in A^-body simulations, such 
as their conditional and unconditional mass functio n, halo 
accreti on histories, merger rates and halo bias (see IZentnerl 
(120071) for a recent review). However excursion set theory 
describes the collapse as spherical, in its o riginal formula- 
tion, o r as ellipsoidal, in the extension due to Shet h & Tormenl 
(1999). This is clearly an important oversimplification of 
the actual complex dynamics and, as a result, while quali- 
tatively the prediction of excursion set theory agree with N- 
body simulations, at the quantitative level there are important 
discrepancies, and dynamical evidence in favor of excursion 
set theory, at least in its present formulation, is quite weak 
dRobertson et al. 2008). A related concern is that numerical 
simulations show that there is not a good correspondence be- 
twe en peaks in the i nitial density field and collapsed halos 
(see lKatz et al.1 d 19931) for an early result). 

In this paper we contin ue the investigation of excur sion set 
theory that we started in Maggior e & Riottol (12009 a!) (here- 
after paper I). In paper I we have shown how excursion set the- 
ory can be put on firmer mathematical grounds, and we have 
been able to take into account analytically the non-markovian 
contribution to the evolution of the smoothed density field, 
due to the use of a tophat filter in coordinate space. In the 
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present paper we turn to a reexamination of the physics be- 
hind excursion set theory and we propose a generalization of 
the theory, based on the idea that the critical value for collapse 
of the smoothed density field should be treated as a stochastic 
variable. We discuss how this stochasticity originates phys- 
ically and we show that supplementing excursion set theory 
with a diffusing barrier allows us to capture at least some of 
the complexity of the actual halo formation process, which is 
lost in the spherical or elliptical collapse model. 

Our notation is as in paper I. Namely, we consider the den- 
sity contrast <5(x) = [p(x)~ p]/ p, where p is the mean mass 
density of the universe and x is the comoving position, and 
we smooth it on some scale R, defining 



S(X,R): 



d 3 x'W(\x-x'\,R)6(x'), 



(1) 



with a filter function W(\x — x'\,R). For gaussian fluctuations, 
the statistical properties of the fundamental density field 6(x) 
are embodied in its power spectrum P(k), defined by 



(m~S(k'))=(2TT) 3 6 D (k + k')P(k), 



(2) 



where 6(k) are the Fourier modes of <5(x). From this one finds 
the variance a 2 (R) of the smoothed density field 



a 2 (R) = (6 2 (R)) . 



(3) 



If we smooth the density field with a tophat filter function in 
coordinate space, the mass M associated to a smoothing radius 
RisM=(4/3)irR 3 p, and we can consider a as a function of 
M, rather than of R. The ambiguities involved in assigning a 
mass M to a smoothing scale R when one uses a different filter 
function have been discussed in detail in paper I. 
The halo mass function dn/dM can be written as 



dn(M) 



p d\na~ l (M) 



(4) 



dM ' M 2 d\nM 

In Press-Schechter theory dPress & Schechterl ll974l) and in 
excursion set theory theory (Bon d et al.lll991l) the function 
f(a) is predicted to be 



/psO): 
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(5) 



where 6 C ~ 1 .686 is the critical value in the spherical collapse 
model. This result can be extended to arbitrary redshift z by 
reabsorbing the evolution of the variance into 6 C , so that 8 C 
in the above result is replaced by S c (z) = 5 C (0)/D(z), where 
D(z) is the linear growth factor. However, eq. Q is valid only 
if the density is smoothed with a sharp filter in momentum 
space, and in this case there is no unambiguous way of as- 
signing a mass to a region of radius R. In paper I we have 
been able to extend this result to a tophat filter in coordinate 
space. In this case the computation is considerably more diffi- 
cult. In fact, when the density perturbation is smoothed with a 
sharp filter in momentum space, 5(R) obeys a Langevin equa- 
tion with respect to the "pseudotime" variable S(R) = o (R), 
with a Dirac delta noise. This means that the dynamics is 
markovian, and that the probability II(<5, 5) that the density 
contrast reaches the value S at "time" S satisfies a Fokker- 
Planck (FP) equation, with an "absorbing barrier" boundary 
condition II(i5 e ,5) = 0. For different filters the dynamics be- 
comes non-markovian, and 11(6, S) no longer satisfies a local 
diffusion equation such as the FP equation. In paper I we 
have been able to formulate the problem of the computation 
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of 11(5,5) in terms of a path integral with boundaries and we 
have found that the result can be split into a "markovian" and 
a "non-markovian" part. The markovian part simply gives 
back eq. ©, where now a 2 (M) is the variance computed with 
the tophat filter in coordinate space, while the non-markovian 
terms can be evaluated perturbatively. To first order, we found 

/(,) = (i- K) p) 1/2 ^-^) + 4= s ^r(o,fi- 2 

\7T J a V27T cr V 2fJ 

where 

k(R)= lim ^fj*^ - 1 ~ 0.4592-0.0031 fl, (7) 

R'^ac {d z (R')) 

R is measured in Mpc//z, T(0,z) is the incomplete Gamma 
function, and the numerical value of k(R) is computed as- 
suming a ACDM model compatible with the WMAP 5yrs 
data and a tophat filter function in coordinate space. Ob- 
serve that for a sharp filter in momentum space, (S(R')S(R)) = 
(6 2 (ma\(R,R')) so k(R), as defined by the first equality in 
eq. (J7j», vanishes. 

However, neither eq. (O nor eq. © perform well when 
compared to cosmological A^-body simulation. Indeed, PS 
theory predicts too many low-mass halos, roughly by a fac- 
tor of two, and too few high-mass halos: at a~ l = 3 (high 
masses correspond to small values of a), PS theory is already 
off by a factor O(10). The mass function given in eq. ©, in 
the interesting mass range, is everywhere lower than the PS 
prediction and therefore, while it improves the agreement at 
low masses, it gives an even worse result at high masses, see 
Fig. 9 of paper I. Thus, it is clear that some crucial physical 
ingredient is still missing in the model. This is not surpris- 
ing at all, given the use of the simplified spherical collapse 
model. In the large mass limit the result cannot be improved 
by turning to the ellipsoidal collapse model since, as we will 
review in the next section, at large masses the barrier for ellip- 
soidal collapse reduces to the one for spherical collapse. The 
aim of this paper is to show that treating the collapse barrier 
as a stochastic variable allows us to capture some of the com- 
plicated physics that is missed by the spherical or ellipsoidal 
collapse model, and we will see that this modification gives 
just the required behavior in the large mass limit. 

The organization of the paper is as follows. In Section [2] 
after recalling how a moving barrier em erges from the ellip- 
soidal collapse model (S heth et ai1l2001l) . we discuss in detail 
the physical motivations for the introduction of a stochastic 
barrier. In Section[3]we compute the halo mass function with 
a diffusing barrier, both for the markovian case and including 
the non-markovian corrections, and in Section|4]we compare 
our prediction for the mass function with A^-body simulations. 
Section|5]contains our conclusions. 

2. THE ELLIPSOIDAL COLLAPSE BARRIER AND THE 
DIFFUSING BARRIER 

The fact that extended PS theory gives a qualitatively 
correct answer but fails at the quantitative level has led 
many authors either to resort to fits to the Af-body simula - 
tions, see e.g . ISheth & Tormenl (119991): ISheth et al.1 (120011): 
Ljenkins et all (120011): IWarren et al l (120061): iTinker et al l 
d2008l) : lPillepichetal . (2008): iGrossi et al l (120091), or to look 
for imp rovements of the spherical collapse model. ISheth et all 
(2001) took into account the fact that actua l halos are triaxial 
Bardeen et"aT1ll986t iBond & Mversll 19961) . and that the col- 
lapse of halos occurs along the principal axes. As a result, the 
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ellipsoidal collapse barrier B acquires a er-dependence, 

1.231 



B(a) ~ 8 C 



1+0.47 



a 



(8) 



Physically this reflects the fact that low-mass halos (which 
corresponds to large a) have larger deviations from sphericity 
and significant shear, that opposes collapse. Therefore low- 
mass halos require an higher density to collapse. In contrast, 
very large halos are more and more spherical, so their effec- 
tive barrier reduces to the one for spherical collapse. 

It is apparent that the use of a moving barrier of the form 
©, by itself, cannot improve the agreement with A^-body sim- 
ulations in the large mass limit since, for large masses (which 
correspond to a — > 0), B(S) reduces to the value for the spheri- 
cal collapse and therefore we get back the incorrect prediction 
of extended PS theory. More generally, since the barrier is re- 
ceding away from its initial location 8 C , it is more difficult 
for the smoothed density perturbation to reach it, at any a, so 
the use of eq. ([8]i simply gives a halo mass function which is 
everywhere smaller than the PS prediction. 

In order to improve the agreement between the prediction 
from the excursion set me thod with a n ellipsoidal collapse and 
the TY-body simulations, Shet h et ail (1200 ll) found that it was 
necessary to introduce a new parameter a (which, when they 
require that their mass function fits the GIF simulation, turns 
out to have the value a ~ 0.707, i.e. y/a ~ 0.84) and postulate 
that the form of the barrier is rather 



B(a) ~ \fa 8 C 



1+0.47 



fade 



1.23 



(9) 



It is important to stress that, in Sheth et al. (2001), the param- 
eter a is not derived from the dynamics of the ellipsoidal col- 
lapse. Rather on the contrary, the ellipsoidal collapse model 
predicts a = 1 because in the limit a — > the barrier mus t re- 
duce to that of spherical collapse. In ISheth et alj d2001l) the 
parameter a is just introduced by hand in order to fit the N- 
body simulations. 

To clarify the origin of this parameter, it is useful to recall 
how eq. ([8]l emerges. One consider s the gravitational co llapse 
of a homogeneous ellipsoid, as in lBond & Mversl dl996l) . De- 
noting by <fi the peculiar gravitational potential at the location 
of an ellipsoidal patch, the deformation tensor is djdi4>, and 
its eigenvalues Ai , A2, A3 (ordered so that Ai < A 2 < A3) char- 
acterize the shape of the ellipsoid. In the linear regime, using 
Poisson equation, the density contrast 8 is given by the trace 
of the deformation tensor, so 8 = Aj + A2 + A3. In a gaussian 
random field, the probabil ity distribution of th e eig envalues is 
known , and is given by dDoroshkevichl (fl970); see lLam et al.l 
(2009) for a recent generalization to the non-Gaussian case) 



p .(Ai,A 2 ,A 3 ) = 



15 3 



-(A 2 -A 1 )(A 3 -A 2 )(A3-A 1 ) 



x exp 



37? 15h\ 
a 2 2a 2 J 



(10) 



where /; = Ai + A 2 + A3 = 8 and 7 2 = AiA 2 + A 2 A3 + A1A3. 
By integrating out Ai and A 2 at fixed 8, with the constraint 
Ai < A 2 < A3, one verifies that 8 has a gaussian distribution, 
with variance a 2 . Rather than using the three eigenvalues as 
independent variables one can use 8, together with the ellip- 
ticity e and prolateness p, defined by 

^3~Ai , f n 



Ai + A3 — 2A2 
28 



(12) 



From p a (A\ , A 2 , X^dXxd^dX?, one can derive the distribution 
probability g a (e, p\8)dedp for e.p at fixed 8. The result is 
( Bardeenetklll986t ISheth et al.ll200Tl) 



g< T (e,p\S) = 



exp 



"2 a- 



Oe 2 + p 2 ) 



(13) 

To define a barrier one needs a c riteriu m for collapse in the 
ellipsoidal case. In She th et alj (UpOl) collapse along each 
axis is stopped so that the density contrast at virialization is 
the same as in the spherical collapse model, i.e. 179 times 
the critical density of the universe. Given this criterium, the 
critical value for collapse in the ellipsoidal model, 8 e \\, is a 
function of e,p, w ell approximated by the implicit relation 
(ISheth et alj|200lh 



5 e ii(e,p) 



l+(3 



5(e z ±p z ) 



81 



(14) 



where 8 C is the critical value for spherical collapse, the plus 
(minus) sign holds for p negative (positive), (3 ~ 0.47 and 
7 ~ 0.615. The barrier ([8]l follows if one replaces e and p 
with their most probable values according to the distribution 
g(e,p\S), which are e = ct/(<5v / 5) and p = (and furthermore 
one replaces 8 e u(e,p), on the right-hand side of eq. (TBI , with 
S). 

As we already mentioned, this barrier always stays above 
the spherical collapse barrier at 8 C , and reduces to the spher- 
ical collapse one for a = 0, i.e. for large masses. For our 
purpose, it is however important to note that this result only 
holds if e and p are replaced by their most probable values e 
and p. For generic values of e and p the critical value for col- 
lapse can be either higher or lower than 8 C . This results in a 
"fuzzy" threshold d Audit et. alj|1997t iLee & Shand"arinll 19981: 
ISheth et ai]|2001l) . with a probability distribution that extends 
even to values smaller than 8 C . To compute the variance of 
the barrier due to this effect we use a slightly more accurate 
expression for the critical value of the ellipsoidal co llapse as 
a function of the eigenvalues A, dSandvik et al. 20061). 

MAi , A 2 , A 3 ) = 5 e [1 + a,(A 2 - Aj ) Q2 + a 3 (2A 3 -A 2 -Ai)" 4 ] , 

(15) 

where a x ~ 0.2809, a 2 = 1.3557, a 3 ^ 0.070, a 4 ~ 1.41205. 4 
We stre ss that this expression for the critical value of 5 was 
found in lSandvik et al.l (|2006) requi ring an accura t e rep resen- 
tation of the ellipsoidal collapse o flBond & Mversl d 19961) . and 
not by fitting to A^-body mass functions. The average value of 
this barrier is 



Bell(<r) = <£ell(Ai,A2,A 3 )) 
A 3 /-Aj 



(16) 



dXi 



dX~ 



dX\ <5eii(Ai, A2,A3)/? CT (Ai, A 2 , A3), 



where p CT (Ai,A 2 ,A3) is the probability distribution given in 
eq. ( TTOb . In Fig.Q]we compare this expression for B e \\{cr) with 
the expression given in ([8]). We see that they provide two sim- 
ilar representation of the average barrier for ellipsoidal col- 
lapse. The variance of B e n(a) is given by 

4 As in eq. (8), we are considering the barrier at z = 0, when the linear 
growth factor D(z) = 1 . The expression for gen eric z is obtained by r eseating 
<5 r -» S c /D(z) and A, -+ XjD(z), see eq. (18) of lSandvik et"aIU 2006). 



4 



3.0 r 




1.5 - 



1.0 - 
(1.5 - 

0.0 t , , , , 1 , , , , 1 , , , , 1 , , , , 1 

0.0 0.5 1.0 1.5 2.0 

a 

FIG. 1. — Two possible representations for the ellipsoidal collapse barrier 
as a function of <r. The red dashed line is the barrier g iven in eq. (8). The 
blue solid line is the function B c ll("") obtained from eqs. i 1 5 i and \Y6\ . 
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FIG. 2. — The function B c ll(o") (blue solid line), together with the curves 
B e ii(cr) ± 3£g d] (cr) (dashed lines). The horizontal red solid line is the spher- 
ical collapse barrier B = S c . 



Y, 2 B Ja)=([S ell (\ u \2,h)-B ell (<j)] 2 ) (17) 

~c5 2 [0.00805cr 2Q2 + 0.00489cr 2Q4 + 0.00305cr Q2+Q4 ] , 

where again the averages have been computed using the dis- 
tribution ( [Tol l. In Fig. |2 we compare B e n{a) with the curves 
-B e ii(er) ± 3£b c1I (i7). We see that a fluctuation of the barrier at 
the 3Sb l11 level can bring the threshold for collapse well be- 
low the constant value 8 C derived from the spherical collapse 
model (see also Fig. 7 of lSheth et ail d2001l) ). 

This result already makes it clear that, as a matter of princi- 
ple, the critical value for collapse is unavoidably a stochastic 
variable. However, the fluctuations of the barrier discussed 
above by no means exhaust all possible sources of stochas- 
ticity in the actual physical problem. For instance, halos 
are subject to tidal effects due to their environment, which 
also results in a distribution of values for the collapse barrier 
(Desiacques 2008). More generally, modeling dark matter ha- 
los as smooth and homogeneous ellipsoids characterized by 
the eigenvalues A,-, even when taking into account their dis- 
tribution probability, is still a significant oversimplification. 
For instance, a patch that is collapsing might have significant 
non-linear substructures, whose presence influences its criti- 
cal value for collapse. All these effects contribute to the scat- 
ter of the values of the threshold for collapse. 

Last but not least, the very definition of what is a dark 
matter halo is a non-trivial problem both in numerical sim- 
ulations and observ ationally (for cluster observations, see 
iJeltema et al.1 (120051) and references therein). In simula- 
tions, halos are usually identified either through a Friends-of- 
Friends (FOF) algorithm, or using spherical overdensity (SO) 
finders. However actual halos are triaxial, rather than spher- 



ical, and often messier than that, and there is nothing funda- 
mental or rigorous in either choice, both being largely a matter 
of convenience. FOF halo finders track isodensity profiles and 
might be more relevant for Sunayev-Zeldovich or weak lens- 
ing, while SO finders may be more relevant for cluster work. 
Searching for halos using for instance a spherical overdensity 
finder, when halos are at best triaxial and often more irregular, 
introduces a further source of statistical fluctuations, both in 
the number count of halo, and in the assignment of the mass. 
A similar concern is that the exact definition of a virialized 
halo depends on the what one means exactly by "virialized". 

So in the end, in a given A^-body simulation, each patch of 
the initial density field that eventually collapses to form a halo 
at a given epoch, has a smoothed overdensity that does not 
have in general exactly the value predicted by the ellipsoidal 
collapse model, but rather fluctuates around it with fluctua- 
tions that are determined by various factors, such as the dis- 
tributions of the eigenvalues of the deformation tensor, the 
details of the halo finder algorithm, or other details related 
to the environment, the presence of non-linear substructures, 
etc., as discussed above. 

Motivated by these considerations, we propose in this pa- 
per to extend excursion set theory by considering a first- 
passage time problem in the presence of a barrier that fluctu- 
ates stochastically. Given that the fluctuations in the collapse 
threshold depend, among other things, on the exact details of 
the halo definition (halo finder, virialization critierium, etc.), 
the mass function computed from excursion set theory with 
such a stochastic barrier will depend on these details, too. 
This is a positive aspect because the actual halo mass func- 
tion o btained from jV-body simulations depends on the halo 
finder (White 2001). For instance, with FOF finders the mass 
function depends on the link-length used, and in particular the 
value a ~ 0.707 given in eq. © holds for a link-l ength equal 
to 0.2 times the mean inter-particle separation (Sheth et al. 
l200l . 5 When halos are identified with a SO finder, the mass 
functio n depends on t he val ue chosen for the overdensity A, 
see e.g. iTinker et al.l (120081) . These effects cannot be repro- 
duced by excursion set theory with a barrier which is fixed 
uniquely by the dynamics of the spherical or ellipsoidal col- 
lapse, and which therefore is insensitive to these details. In 
this paper we explore whether the use of a stochastic barrier 
allows us to incorporate into excursion set theory, at least at 
the level of an effective description, a part of the stochasticity 
intrinsic to the actual physical problem of halo formation, and 
due both to the complicated underlying dynamics and to the 
choices that one has made when giving an operative definiton 
of dark matter halos. 

Ideally, one would like to compute theoretically the fluctu- 
ation properties of the barrier. For some effect, such as that 
due to the distribution of eigenvalues of the deformation ten- 
sor, this is possible, as we saw in eq. (1171 1. Unfortunately, 
other effects such as the scatter of the barrier due to the de- 
tails of the halo finder (which, as it will turn out, give a con- 
tribution that dominates over that in eq. ( fl7l i) are much more 
difficult to predict theoretically. We will therefore take a more 
phenomenological approach. We will consider a barrier that 
performs a random walk with a diffusion coefficient Db- At 
least at a first level of description, all our ignorance of the 

5 Observe that, for this link length, a sizable fraction of halos have major 
non-spherical substructures and a signi ficant con tribution to the halo mass 
arises from outside the "virial" radius (Lukic et al. 2009). This underlines 
again the importance of the details of the definition of what is a halo. 
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details of the dynamics of halo formation is buried into this 
coefficient. Solving the first-passage time problem with such 
a barrier we will find that the net effect is that in the halo 
mass function predicted by Press-Schecther or excursion set 
theory the exponential factor changes from exp{-6 2 /(2a 2 )} 
to exp{-a<5 2 /(2er 2 )} (and more generally we must replace ev- 
erywhere S c — > S c y/a), where 

1 



(18) 

see Sec tion [3] This is jus t the replacement t hat w as postu- 
lated in ISheth et al.l (120011) : ISheth & Tormenl (1 19991) in order 
to fit the data. We therefore discover that the Sheth and Tor- 
men (ST) mass function (at least in the large mass limit), is 
just the mass function obtained by excursion set theory with a 
diffusing barrier. 

Having obtained a physical understanding of the parameter 
a that appears in the ST mass function, one can ask whether 
it is p o ssible to go beyon d the approach in Sheth & Tormen 
( 1999): ISheth et al.ld2001l) . where a is simply treated as a fit- 
ting parameter, and try to predict it, by computing the diffu- 
sion coefficient Db- Given that a first-principle computation 
of Db seems difficult, in Section |4] we will turn to iV-body 
simulations themselves. We will see that recent numerical 
studies of the properties of the collapse threshold allows us to 
deduce the value of Db- Given this input we will then com- 
pare our prediction (TT~8T > to the slope of the exponential fall 
of of the mass function, and more generally (including in the 
halo mass function also the non-markovian corrections com- 
puted in paper I) we will compare our analytic form for the 
mass function to the numerical data. Even if in this way we 
must use an input from A^-body simulations themselves, still 
the comparison is quite non-trivial, since the relation ([T8T > is a 
specific prediction of our model. 

3. THE HALO MASS FUNCTION IN THE PRESENCE OF THE 
DIFFUSING BARRIER 

In order to illustrate the idea in a simple mathematical set- 
ting, we consider a barrier that fluctuates over the constant 
spherical collapse barrier B = 8 C . In the large mass limit the 
ellipsoidal collapse barrier reduces to the spherical collapse 
barrier, so we expect that this approximation should be ade- 
quate for computing the effect of a stochastic barrier on the 
high-mass tail of the halo mass function. We also consider 
a barrier B that fluctuates in such a way that its mean root 
square fluctuation depends linearly on the variance of the 
smoothed density field <r{R), 

T, B (a(R)) = ((B-(B)) 2 y/ 2 = a(R)VDB~ 1 (19) 
where D B is a numerical coefficient. This choice is partly mo- 
tivated by mathematical simplicity. Furthermore, we will see 
in SectionHthat there is some evidence from Af -body simula- 
tions that for small <r(R) the barrier diffuses just as in eq. ( fT9l . 

This form of the barrier corresponds to a brownian motion. 
Recall in fact that, if a particle performs a one-dimensional 
brownian motion, its position x(t) has a variance given by 
((x(t)-xo) 2 } 1 ' 2 = \Dt, where D is the diffusion coefficient. 
In the excursion set method S(R) = cr 2 (R) plays the role of 
a "pseudotime" variable, so eq. (TT9b means that the collapse 
barrier performs a brownian motion around its initial posi- 
tion (B), with a diffusion coefficient Db- We will refer to 
the model in which the barrier's scatter behaves as in eq. ( fT9l 
as a "diffusing barrier". As we will see below, the first pas- 
sage time problem in the presence of a diffusing barrier can be 



solved analytically, so eq. (T% provides at least a useful toy 
model for understanding the effect of a more general stochas- 
tic barrier. We emphasize however that the idea that we are 
proposing is more general, and can in principle we imple- 
mented for a generic functional form of the barrier B(a) and 
of its fluctuation Sb(ct), although the associated first-passage 
time problem becomes more complicated, and might require 
the numerical generation of an ensemble of trajec tories using 
Monte Carlo simulations, as in BonditaD Jl99lh . 

Intuitively, we can understand why a diffusing barrier can 
help to reproduce the numerical A^-body results. We are in- 
terested in events, corresponding to cluster masses, that arise 
from rare fluctuations, on the far tails of the probability dis- 
tribution. For instance at er" 1 = 3, the PS theory prediction for 
f(a) is about 10~ 5 , and we are searching for a mechanism that 
brings this number up to the observed value of about 10~ 4 . 
Even if, on average, a barrier has equal probabilities of fluc- 
tuating toward values lower that 6 C as toward values higher 
than 6 C , still the fluctuations of the barrier toward lower values 
can have a much more significant effect (consider for instance 
what happens to a dam if on a rare occasion it is lowered). In 
fact, this is true even if most of the flucutation where above 
6 C , and only rare flucutations were below 8 C which, as we will 
see, is the case when we consider fluctuations over the ellip- 
soidal collapse barrier, see also Fig. [2] In the analogy with the 
dam, rare occasional lowerings of the dam can produce sub- 
stantial flooding, even if many more flucutations rather raise 
it. 

To verify formally this intuition, we neglect at first the non- 
markovian corrections due to the filter, discussed in paper I. 
We denote by U(B(0),B;S(0),S;S) the joint probability that, 
at "time" S, the barrier has reached by diffusion the value B, 
starting from the initial value B(Q) = 6 C , while the density con- 
trast has reached the value 6, starting from the initial value 
6(0) = 0. In the markovian case the probability distribution 
obeys a Fokker-Planck equation. The fact that the "particle" 
described by 6(S) and the barrier B(S) both diffuse indepen- 
dently means that the joint probability distribution satisfies the 
two-dimensional FP equation 



on _D d 2 n D B d 2 n 
~dS~2~d6 2+ ^'dB 2 



(20) 



In our case the diffusion coefficient of D = 1, see e.g. eq. (20) 
of paper I, while Db is the diffusion coefficient of the barrier. 
To solve this equation it is convenient to introduce a "time" 
variable t = S/5 2 , and the variables 



6r~B 



Xl = 



x z = 




so eq. (f20b becomes 



dU 



1 s 2 n 1 <9 2 n 

■ + ■ 



2 Ox 2 



2 dx 2 



(21) 
(22) 

(23) 



In term of these variables the barrier starts at xi(0) = while 
the "particle" starts at X2(0) = 1. The boundary condition 
is that IL(B(0),B;6(0),6;S) vanishes when 6(S) = B(S), i.e. 
n(xi(0),xi;jC2(0),X2;f) vanishes when \/DbX\ =Xz- We de- 
fine 6 from \[D~b~ = tan 6, so we have a two-dimensional FP 
equation with the boundary condition that II vanishes on the 
line xi = X2Cot#, see Fig. [3] This problem can be solved by 
the method of images (Redner 2001), and the result is given 



FIG. 3. — Mapping of the scaled density perturbation and collapse barrier 
coordinates in one dimension to the plane (x\,X2). The initial position is in 
(xi = 0,xj = 1) (black dot) and its image point is in (xi = sin20.,V2 = —cos 28) 
(white dot). 



FIG. 4.— The function ns m (<5, 5) with D B = 0.25 (blue solid line) compared 
to the standard excursion set result, i.e. n gm (<5, 5) with Dg = (violet dashed 
line), as functions of <5, for S = 1 . 



by a gaussian centered on (x\ = \,x<i = 0) minus a gaussian 
centered on the image point (x\ = sin2f9,x 2 = -cos2f9), 



l 

~ 27Tf 



0,xi;x 2 (0) : 

[* 2 +fe-l) 2 ]/2f . 



: \,x 2 ;t) 

-[Ui-sin2e) 2 +(jE 2 +cos26>) 2 ]/2r 



(24) 



where, as in paper I, we added to II the superscript "gm" to 
remind that this is the solution for gaussian fluctuation with a 
markovian dependence on the smoothing scale 

The probability density for the scaled density perturbation 
to be at the position x% is the integral of the two-dimensional 
density over the accessible range of the scaled collapse barrier 
coordinate x\ , 




A'2 COt 8 



e -fe-«-A Erfc 



dx x W m (x x \x2\t) 

X2 COtf? X 



>2t 



. e - te+ cos20f/2r Erfc ( sin20-x 2 cot0 



lit 



(25) 



where Erfc(z) is the complementary error function and the 
initial conditions xi(0) = and x 2 (0) = 1 are understood. 
Restoring the original variables S, S(S) and B(S), and using 
U(5o;S;S)dS = Il(x2(Oy,X2;t)\dx2\ where \dx2\ = dS/S c , we get 



n gm (5,S) = 



1 



2V2ttS 

_ e -(2,5 £ .cos 2 e-5) 2 /(2S) Erfc 



e " 52 / (2s) Erfc 



cot 9 



(5c- 5) 



5 c sm29-(8 c -5)co\.9 



2S 



,(26) 



where the initial condition So = is understood. The limit of 
non-diffusing barrier is D B — ► + , so 9 — ► + and cot 9 — > +oo. 
Recalling that Erfc(z) — > 2 as z — > -co, we see that in the limit 
Db — * + we recover th e standard r e sult o f excursion set theory 
with a static barrier of iBond et"aT1 (Il991h . 

In Fig. |4]we compare this function, for a diffusion coef- 
ficient Db = 0.25, with the static barrier case. Observe that, 
when Db = 0, the distribution function vanishes for 6 > 5 C , 6 
while for finite Db it is non-zero for all values of 5. Of course, 

6 This only holds because, for the markovian term, we can work directly 
in the continuum limit. As we discussed in paper I, if we compute II sum- 
ming over trajectories defined discretizing time in steps e, there are finite e 
corrections and II(<5.S) non longer vanishes for 8 > S c , even for the static 
barrier. 



this reflects the fact that the barrier can in principle diffuse to 
arbitrarily large values of 6. 
The markovian contribution to the first crossing rate is 



F sm (S) = - 



dS 



OS ■ 



(27) 



The evaluation of this expression can be simplified observ- 
ing that d/dS, when acting on (2ttS)~ 1 ^ 2 exp{-<5 2 / (25)}, is the 
same as (l/2)d 2 /dS 2 , and integrating twice by parts d 2 /dS 2 . 
We then find 

S 2 



VW+Db)S 3 / 2 



exp 



2(1 +D B )S 



(28) 



The function f(<f) is obtained from the first crossing rate using 
f(a) = 2a 2 T(a 2 ), see e.g. Section 2 of paper I, so we get 



1/2 



\fa5, 



£ -aS 2 c /(2a 2 ) 



where 



l+D h 



(29) 



(30) 



This is the crucial result of this section. We see that the 
effect of the diffusing barrier is that the exponential factor 
in the halo mass function changes from e\p{-5 2 /(2a 2 )} to 
exp{-fl(5 2 /(2cr 2 )}, with a given by eq. d30b . and more gen- 
erally everywhere in the mass function S c — > S c y/a. This 
is exactly the m odification which was postulated ad hoc in 
IShethetal.l(l200l . 

In fact, even if the expression for n gm (<5, S) given in eq. d26j > 
is interesting by itself, the result for the first-crossing rate 
could have been obtained directly, without even computing 
explicitly II gm ((5, S), simply by observing that the problem in- 
volving a barrier with coordinate xi and diffusing with a dif- 
fusion coefficient D\, and a particle with coordinate X2, dif- 
fusing with a diffusion coefficient D 2 , can be mapped into a 
one-degree of freedom problem, introducing the relative coor- 
dinate x = X2 — Xi . The resulting stochastic motion is go verned 
by an effective diffusion coefficient D e ff = D\ +D 2 dRednerl 
1200 lb . This point can be easily understood considering a 
Langevin equation for the barrier coordinate x\, 



with 



x\ = T)l(t), 

(vi(t)m(t')) = Dt6(t-t'), 



(31) 
(32) 
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and a Langevin equation for the particle coordinate %i = rfyit) 
with (r]2{t) r l2(t 1 )) = Ih.S(f — t'). Then the relative coordinate 
x = X2~xi satisfies x= 77(f) with 77(f) = 772(f) — 771(f) and, if 771(f) 
and 772(f) are uncorrelated, 



(v(t)v(t')) = MOmC)) + {V2(t)vi(t')) 

= (D l+ D 2 )5(t-t'), 



(33) 



showing that the relative coordinate diffuses with an effective 
diffusion coefficient D^+D^. In our case D\ =D B and£>2 = 1. 

We have repeated the above analysis including the non- 
markovian corrections due to a tophat filter in coordinate 
space, to first order, using the results obtained in paper I. The 
explicit computation is performed in Appendix A. The result 
is 



(34) 

where 

R = TTd b - (35) 

This is our result for the halo mass function. We hasten to 
add that this result only holds in the large mass limit, other- 
wise we must consider fluctuations over the ellipsoidal bar- 
rier, rather than over the constant spherical constant barrier, 
and furthermore we have considered a specific form of the 
barrier variance, corresponding to a random walk. We now 
turn to a comparison of our model with A^-body simulations. 

4. COMPARISON WITH A-BODY SIMULATIONS 

The variance of the collapse barrier in Af-body si mula- 
tio ns has been recentl y studied in Robe rtson et al.l d2008f) (see 
als dDalal etaT] d2008h ). For each halo identified in their TV- 
body simulations at z = 0, they calculated the center-of-mass 
of the halo particles at their positions in the linear field at 
z ~ 10 2 and used the density field, smoothed with a tophat 
filter in real space, to compute the overdensity within a given 
lagrangian radius R. This overdensity is then linearly extrapo- 
lated to z = 0. They find that the distribution of such smoothed 
linear overdensities B(a), at fixed a, is approximately log- 
normal in shape with a width A B given by 



1/2 



R5r 



Ab — 0.3cr. 
In a log-normal distribution one has 



{{B-{B)f)^=(e 



1/2 



(36) 



(37) 



and in our case, for small a, (B) ~ S c ~ 1.68. Observe that 
eq. d37"l i is consistent with our estimate d!71 l. In fact, eq. dTTb 
only includes the scatter due to the fluctuations of the eigen- 
values, so it is a lower bound on the actual scatter of the bar- 
rier that, as we discussed in Section[2] can in principle receive 
contributions from many other effects. As we see from Fig. [5] 
the variance given by eq. d37l i is indeed always above that 
given by eq. dT7l >. 

In a ACDM model, a(M) is such that, for values of M 
corresponding to cluster of galaxies, Q3cr(M) <C 1. For in- 
stance, a(M) = 1 for M ~ lO 14 M /r\ while a( M) = 0.6 
for M c± 10 15 Mq/2 _1 (see e.g. Fig. 1 of the review IZentnerl 
d2007l) ). Therefore in the high-mass range As is small and 
we can expand eq. d37l ), obtaining 



FIG. 5. — The varianc e of the barrier Eg from eq. {37) (blue solid line) 
compared to the estimate {T7} (violet dashed line). 



The inclusion of an overall drift of the barrier, such as that in 
eq. ©, as well as higher order terms in the expansion of the 
exponential in eq. (|37| |. provides terms of higher order in a, 
which are subleading in the large-mass regime. 

Equation (l38l has exactly the form of the diffusing barrier 
given in eq. dl9> . with a diffusion coefficient 



Db — (0.3 <5 C ) ~ 0.25 . 



(39) 



In this case our model predicts a relation, given by eq. d30b , 
between the diffusion coefficient Dg of the barrier and the 
slope of the exponential of the halo mass function in large 
mass limit, i.e. we predict 



1 



-0.80. 



(40) 



T he value ([39b has been deduced from the 7Y-body simulations 
of iRobertson et al.l ((2008), where halos where identified with 
a A = 200 spherical overdensity algorithm. We therefore must 
compare our prediction with the value of a obt ained under 
the sam e conditions. This can be obtained from Tink er et"aT] 
(2008), where the same numerical simulation was used to 
study the halo mass function (and its deviations from univer- 
sality, see below). The authors fit their result with a fitting 
function 



f(o)=A T 



+ 1 



ct/o 



(41) 



(B)A B ~035 c a. 



(38) 



and, for a spherical overdensity A = 200, they findAr = 0.186, 
aj = 1.47, bj = 2.57 and ct — 1.19. We add a subscript T, 
which stands for Tinker et al., to distinguish for instance their 
parameter a T from our parameter a. A first indication of the 
agreement of our prediction with the above fitting formula can 
be obtained by comparing the respective exponential cutoff. 
In terms of their parameter ct, our parameter a is given by 
the combination 2ct/S^. Using their value cj ~ 1.19, one 
has 2ct/S^ ~ 0.837, in good agreement with the value d40b . 
This agreement is a non-trivial result. It is true that, to get 
eq. (|40l i, we used an input from the same 7Y-body simulation, 
namely the scatter of the values of the threshold for collapse, 
from which we deduced the diffusion coefficient Db of the 
barrier. However, given this input our model makes a non- 
trivial prediction for the the numerical value of the parameter 
a that appears in the halo mass function. This is very different 
from fitting a directly to the 7Y-body simulations. In principle, 
the prediction a = 1/(1 +Db) could have given rise to a value 
of a very different from the one extracted directly from 7A- 
body simulations, and this would have falsified the diffusing 
barrier model. 

Of course, given that the functional forms of f(a) in 
eqs. d34l ) and ( HTT i are different (which also implies that the 
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FIG. 6.— Our prediction for f(a) given in eq. {34} ( blue solid line) com- N-body simulations given by eq. {41} (blue solid line), setting re = 0. The 

pared to the fit to N-body simulations given by eq. 441 1 (violet dashed line), dashed line marks the line R=l, 

in a log-log scale. 
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FIG. 7. — The ratio R betwee n ou r prediction {34} for f(cr) and the fit to 
/V-body simulations given by eq. 14 It (blue solid line). The dashed line marks 
the line,R = 1. 



normalization constant A 7- in eq. (1411 is not exactly the same 
as the overall factor that we have in front of the exponential in 
eq. (l34l>). the proper way of performing an accurate compar- 
ison is not in terms of the location of the exponential cutoff, 
but rather directly in terms of the full functions f(a). In Fig. [6] 
we compare our prediction for /(er), given by eq. (f34T >. to the 
function fiTT i representing the fit to the TV-body simulation. 
Observe that the vertical axis ranges over more than three or- 
ders of magnitudes. 

To make a more detailed comparison in Fig. [7] we plot, on 
a linear-linear scale, the ratio R between our prediction for 
f(a) and the Tinker et al. fit to TV-body simulations given 
in eq. (RTt . We see that for all values of <j~ l > 0.3 the dis- 
crepancy between our analytic result and the fit to the TV-body 
simulation is smaller than 20%, and for cr" 1 > 1 it is smaller 
than 10%. Considering that our result comes from an ana- 
lytic model of halo formation with no tunable parameter (the 
parameter a is fixed once Dg is given, and we do not have 
the right to tune it), while eq. fiTJ is simply a fit to the data 
with four free parameters, we think that this result is quite 
encouraging. The numerical accuracy is actually the best 
that one could have hoped for, considering for instance that 
we have neglected second-order non-markovian corrections. 
From eq. (l34l we see that, in the computation of the non- 
markovian effect due to the tophat filter in coordinate space 
in the presence of a diffusing barrier, the actual expansion pa- 
rameter is ii = k/(\+Db) which, using eq. (jTj with/?= lOMpc 
and Db — 0.25, has a numerical value R ~ 0.34. Therefore the 
second-order non-markovian corrections, which are propor- 
tional to k 2 , are expected to be of order 10%. Furthermore, 
as we move toward lower masses the effect of the ellipsoidal 
barrier must become important, while eq. ( [34-b has been ob- 



tained using the spherical collapse model, and the variance of 
the barrier shown in Fig. (0 (solid line) has been approx- 
imated by a straight line in eq. d38l . which again is only valid 
at small a. 

In Fig. [8] we show the result that one obtains for the ratio 
R if one includes the diffusing barrier but neglects the non- 
markovian corrections due to the tophat filter in coordinate 
space, i.e. if one sets n = in eq. ([34). We see that the 
agreement degrades, and the discrepancy becomes of order 
30-50%. Thus, while the largest part of the improvement, 
compared to PS theory, comes from the introduction of a dif- 
fusing barrier (recall that PS theory, which predicts a = 1, is 
off by one order of magnitude in the high-mass limit, see 
e.g. Fig. 1 of paper I), still for an accurate computation it 
is important to include the non-markovian corrections due to 
the tophat filter in coordinate space. Observe also that, in 
the large mass limit, the term proportional to the incomplete 
Gamma function in eq. (l34l > is subleading, and the effect of 
the filter is basically to reduce the the overall numerical fac- 
tor, compared to PS theory, by a factor 1 - k. Note that in 
the ST mass function the numerical value of the overall con- 
stant is fixed by hand, by imposing the normalization con- 
dition that all the mass ends up in virialized objects. In our 
case, in contrast, the mass function comes out automatically 
with the correct normalization, as we already showed in Sec- 
tion 5.4 of paper I. The derivation given in eqs. (126)-(128) 
of paper I goes through trivially when 5 C — » 8 c -\fa, so the term 
proportional to the incomplete Gamma function in eq. d34t 
ensures that the mass function is properly normalized, when 
the amplitude of the term proportional to exp{-aS 2 / (2<r 2 )} is 
reduced by a factor 1 — K. 

5. CONCLUSIONS 

In this paper we have proposed a generalization of excur- 
sion set theory, based on the idea that the threshold for col- 
lapse should be treated as a stochastic variable, fluctuating 
around the ellipsoidal collapse barrier (or, in the large mass 
limit, around the spherical collapse barrier). We have seen 
that fluctuations in the threshold arise naturally from a num- 
ber of physical effects. For instance, even within the highly 
simplified description in which a halo is modeled as a smooth 
and homogeneous ellipsoid, fluctuations in the collapse bar- 
rier arise from the fact that the eigenvalues of its deforma- 
tion tensor are stochastic variables, governed by a distribution 
probability. Only when one averages over this probability dis- 
tribution one recovers a barrier which is a fixed function of 
the variance <r 2 (R) of t he smoothed densit y field. Otherwise, 
as already discussed in lAudit et. al.l dl997l) ; ILee & Shandarinl 
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dl998l) : ISheth et alj (1200 ll) one has a "fuzzy barrier" which 
fluctuates around the ellipsoidal collapse value, with fluctua- 
tions that can occasionally bring the critical value for collapse 
even below the spherical collapse value <5 C , see e.g. Fig. [2] As 
we have discussed, many other effects, such as the details of 
the halo finder, the environment, or the presence of non-linear 
substructures, contribute to these fluctuations. 

For mathematical simplicity in this paper we have restricted 
ourselves to a barrier that performs a diffusive motion, with 
diffusion constant Dg, around the constant value given by 
spherical collapse. We expect this to be a good approxima- 
tion in the large mass limit. For such a barrier we have found 
that the first-passage time problem can be elegantly solved, 
and leads to a very simple result. Namely, in the mass func- 
tion one must replace S c — * y/aS c , where a = 1/(1 + Db). The 
replaceme nt 8 C — > \fa8 c is ju st the modification that was pos- 
tulated in lSheth et ail (120011) . in order to fit the results of N- 
body simulations. The diffusing barrier model therefore offers 
a physical understanding of this modification of the PS mass 
function. 

We have then combined our diffusing barrier model with 
the non-markovian corrections due to a tophat filter in coor- 
dinate space computed in paper I, and we have presented an 
analytic expression for the halo mass function, valid for large 
masses. This result can be compared with existing A^-body 
simulations. We have inferred the va lue of Db from the results 
presented in lRobertson et al.1 (12008). and given Db our model 
predicts the corresponding value of a for the same simula- 
tion. Our mass function, with a fixed in this way, is then com- 
pared to the corresponding A^-body simulations in Figs. [6]T7] 
The agreement is better than 20% for all er" 1 > 0.3 (corre- 
sponding approximately to halo masses M> lO n M //i) and 
better than 10% for all a~ l in the interval 1 < a~ l < 3 (corre- 
sponding approximately to halo masses from M ~ 1O 14 M / h 
toM~ 10 I5 M Q //i). 

We conclude with an assessment of what can be obtained 
from excursion set theory, when it is combined with a diffus- 
ing barrier model, and with a discussion of possible improve- 
ments of the model. First of all one should stress that this 
theoretical model, using relatively simple ingredients, investi- 
gates a very complex phenomenon such as halo formation. It 
is therefore encouraging that it nevertheless provides an ana- 
lytic result for the halo mass function that agrees with the N- 
body data with a precision better than 20% over four decades 
of halo masses (and the precision becomes of order 5-10% in 
the higher mass range). Considering that over this mass range 
the halo mass function changes by more than three orders of 
magnitude, this is a non-trivial result. 

We also stress that, when comparing our result to the N- 
body data as in Fig. [6] we had no freedom of adjusting free 
parameters. The functional form of the halo mass function 
was derived from our model, and in this sense it has a dif- 
ferent meaning compared to many fitting functions that have 
been proposed in the literature with the only aim of repro- 
ducing the A^-body data. We needed an input from Af-body 
simulations, namely the scatter of the values of the threshold 
for collapse, from which we deduced the diffusion coefficient 
Db of the barrier. However, given this input our model makes 
a prediction for the the numerical value of the parameter a 
that appears in our halo mass function. This is different from 
fitting a directly to the A^-body simulations. The fact that the 
resulting halo mass function agrees with A^-body simulations 
much better than the original extended PS theory lends sup- 



port to the idea that the diffusing barrier model provides an 
effective way of including, within the excursion set theory 
framework, a number of physical effects that are lost when 
excursion set theory is combined with the simpler spherical 
or ellipsoidal collapse models. 

For precision cosmology, especially in the future, an accu- 
racy such as the one that we have achieved is probably not yet 
sufficient. Of course, generally speaking, analytic models are 
not meant to compete with very time-consuming numerical 
simulations as far as accuracy is concerned. Rather, their role 
is to provide some physical understanding and some guidance. 

However, improvements of the model are certainly possi- 
ble. In particular, rather than considering a diffusive motion 
around the constant value S c , one should consider the actual 
behav ior of the thresho l d and of its variance with a. Accord- 
ing to iRobertson et al.1 (120081) . the average value of the bar- 
rier basically follows the prediction dHJ of the ellipsoidal col- 
lapse, and the scatter around it has a log-normal distribution. 
In general, fluctuations of the barr ier below 5 C are rare ( the 
vast majority of points in Fig. 3 of Robertso n et al.l (120081) lie 
above S c , since the average value follows the ellipsoidal col- 
lapse barrier, which is a rising function of a). However, as we 
discussed in Section [3] one should not forget that the fluctu- 
ations leading to massive clusters are rare events, which be- 
long to the high-mass tail of the distribution function. Since 
the collapse of a halo depends exponentially on the square 
of the height of the barrier, even the rare occasional fluctu- 
ations that bring the barrier below S c can end up enhancing 
significantly the formation probability of the rarest objects. 
The first-passage time problem with such a stochastic barrier 
might be hard to solve analytically, but one could simply inte- 
grate nume r icall y the corresponding Langevin equation, as in 
iBond et all (1199 ll) . In the limit of small a one should recover 
the results presented in this paper, but for intermediate values 
of a there will be corrections. We plan to investigate this issue 
in future work. 

Another possible future developement is the investiga- 
tion of whether universality violations can be accounted for 
within the excursion set theory framework, supplemented by 
a stochastic barrier. Recall that, within excursion set theory, 
the mass function can be written as in eq. (|4J), where the func- 
tion / depends only on the variance a of the smoothed den- 
sity field. Thus, this function has a universal form in the sense 
that its dependence on redshift and on cosmology enters only 
through the dependence of the variance a(R,z). In A^-body 
simulations there are indicatio ns of violations of universal- 
ity, at approximately 10% level dReed et al.l20 06: Tink er et"aT] 
120081) . The evolution with red shift of the exponential cutoff 
is minimal dTinker et al.l 12008b while a redshift dependence 
shows up in the coefficients aj, bj and Aj of eq. (HTI) . Since 
cj and therefore our parameter a, do not show appreciable 
dependence with z, within our model these violations of uni- 
versality cannot be ascribed to a redshift dependence of the 
diffusion coefficient Db- However, Db only reflects the scat- 
ter of the barrier near a = 0. In excursion set theory, the a- 
dependent prefactors in front of the exponential, in the halo 
mass function, rather originate from th e shape of the el lip- 
soidal collapse barrier away from a = (ISheth et aLCOOlh . It 
would be interesting to study, with A^-body simulations, the 

7 Observe however that these violations of universality depends on the halo 
finder algorithm, and at least with some halo finders they can b e accounted 
for by systematic corrections due to the finite simulation volume ILukic et"all 
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shape and the scatter of the collapse barrier as a function of 
the redshift at which the halos eventually collapsed and viri- 
alized, i.e. to repeat fo r different z the analysis performed in 
iRobertson et al.l (120081) at z = 0. One could then generate an 
ensemble of trajectories by integrating numerically the cor- 
responding Langevin equation, and study when these trajec- 
tories first pierce such a stochastic barrier. In this way one 
could investigate whether excursion set theory supplemented 
by a stochastic barrier can account for the observed deviations 
from universality. The interest of such a procedure is that, 
if this were indeed the case, one would have obtained some 
insight into the physical mechanisms responsible for the vio- 
lations of universality. If, in contrast, this procedure should 
not reproduce the observed universality deviation, one would 
have to conclude that this is an intrinsic limitation of excur- 
sion set theory. In any case, one would get a better under- 
standing of what can, and what cannot, be explained within 
the framework of this theoretical model. 



Finally, another interesting test of our model that could be 
performed with A^-body simulations is the computation of the 
barrier scatter, and hence of Db, with different halo finders 
(i.e. with different values of A in the spherical overdensity al- 
gorithm, and with different link-length in the FOF algorithm). 
Changing the halo finder changes the exponential factor in 
the mass function, i.e. the constant a, and the diffusing barrier 
model predicts that the diffusion coefficient Db should change 
accordingly, in such a way that the relation a = 1/(1 +Db) is 
preserved. 

We thank Sabino Matarrese, Ravi Sheth, Cristiano Porciani 
and Sidney Redner for useful discussions. The work of MM 
is supported by the Fond National Suisse. The work of AR is 
supported by the European Community's Research Training 
Networks under contract MRTN-CT-2006-035505. 



APPENDIX 

A. INCLUSION OF NON-MARKOVIAN CORRECTIONS 

In this appendix we compute the effect of the non-markovian corrections due to a tophat filter in coordinate space. As it was 
shown in paper I, when we use a tophat filter in coordinate space the two-point correlation function can be written as 

{5(S0S(S2)} = min(S,, S 2 ) + A(S h S 2 ), (Al) 

where 

A(S u S 2 )~k - , (A2) 

and k ~ 0.45. The first term in the right-hand side of eq. ( IA1 b is responsible for the markovian contribution to the dynamics, and 
it originates from a Dirac-delta gaussian noise; the second term provides the non-markovian contribution. The reader is referred 
to paper I for more details. 
The fact that the barrier diffuses with a diffusion coefficient Db means that 

(B(S l )B(S 2 )) = D B min(5, , S 2 ) • (A3) 

More generally, even the motion of the barrier can be subject to non-markovian effects, so eq. (IA3b should be generalized to 

{BiSOBiSi)) = D B mm(S u S 2 ) + A B (S U S 2 ) ■ (A4) 

Making the rather natural assumption that (5(Si)B(S 2 )) = and introducing the variable X(S) = 5(S)-B(S), we see that 

(X(S l )X(S 2 )) = (l+D B )mm(S u S 2 ) + A(S u S 2 ) + A B (S u S 2 ). (A5) 

Thus, our problem becomes formally identical to a problem for a single degrees of freedom X(S), with an absorbing boundary 
condition at X = 0, with diffusion coefficient (1 +D B ), and non-markovianities described by A(S\ ,S 2 ) + A B (Si,S 2 ). 

We now make the assumption that Ab{Si,S 2 ) is small with respect to A(S\,S 2 ). This assumption could be tested by ex- 
tracting the correlato r (B(S\)B(S 2 )) from the A^-body simulations, similarly to how the variance (B 2 (S)) has been computed in 
IRobertson et al.1 (120081) . The effect of a non-vanishing As can be included perturbatively using the technique that we developed 
in paper I, just as we did for A{S\,S 2 ). (Actually, we expect that the two-point function of the critical value for collapse B{S) 
receives non-markovian corrections due to the same smoothing procedure that gives the non-markovian corrections to S(S) so, 
if this is the dominant effect, a plausible expectation is that (B(Si)B(S 2 )) = D B [mm(Si,S 2 ) + A(S\,S 2 )], i.e. the barrier has the 
same two-point function as S{S), apart from the overall diffusion coefficient, so A B (S\ , S 2 ) = D B A{Si , S 2 ). If this is the case, R in 
eq. (f35T > below is replaced by n. This would entail a 0(25)% modification of the non-markovian correction computed below.) 

When Ab can be neglected, the computation of the halo mass function to first order in the non-markovian corrections can be 
performed introducing a rescaled "time" variable S = (1 +Db)S. Then, using the explicit expression (IA21 >. we get 

(X(S l )X(S 2 ))=mm(S u S 2 ) + k-^-l - (A6) 

S 2 

where R = This is the same problem that we have already solved in paper I, with n replaced by R and S replaced by S, so the 
solution can be written immediately, and is given by eqs. fl34"l > and ( 1301 ). 
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